#########################################################################   
#                                 INFO                                  #   
#########################################################################  

# PROJECT: Exposure to State Violence and Substance Use
# AUTHOR: Alexandra Blackman; Elizabeth Nugent; Sarah Kammourh
# PURPOSE: analysis of drug attitudes
# CREATED: October 2022
# INPUTS:  egy_qual_vars.xlsx

#########################################################################   
#                                 SETUP                                 #   
######################################################################### 

######## ENVIRONMENT

rm(list = ls()) 
setwd("~/Dropbox/Egypt_Alex Sarah Liz/")

######## PACKAGES
need <- c("doBy","foreign", "lubridate","extrafont", "car", "plyr",
          "reshape", "reshape2", "openxlsx", "tidyr", "dplyr", "ggplot2", "ggthemes",
          "RColorBrewer", "rms", "lmtest", "sandwich", "msm","stargazer","texreg",
          "readr","purrr","tibble","stringr","forcats","jsonlite","gridExtra",
          "xml2","httr","rvest","DBI","hms","blob","magrittr","glue", "stringr",
          "readxl","tidyverse","multiwayvcov","readstata13","ClustOfVar","scatterplot3d",
          "mlogit","MASS") # packages needed
have <- need %in% rownames(installed.packages()) # installed packages
if(any(!have)) install.packages(need[!have]) # install missing packages
invisible(lapply(need, library, character.only=T)) # load needed packages


######## SET ALL TABLES TO SHOW "NA" BY DEFAULT

table = function (..., useNA = 'ifany') base::table(..., useNA = useNA)  

######## CLUSTERING FUNCTION

super.cluster.fun <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}

f.test <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  ftest <- waldtest(model, vcov = vcovCL, test = "F")
  return(ftest)
} 


#########################################################################   
#                 LOAD DATA                                             #   
######################################################################### 

# Load survey data
egy_qual <- read.xlsx("~/Dropbox/Egypt_Alex Sarah Liz/Data/egy_qual_vars.xlsx")
head(egy_qual)

#########################################################################   
#                   EXPOSURE TO DRUG USE                                #   
######################################################################### 

# Tramadol Usage
table(egy_qual$tramadol)
egy_qual$tramadol[egy_qual$tramadol == ""] <- NA
egy_qual$used_tramadol <- ifelse(egy_qual$tramadol=="نعم",1,0)
table(egy_qual$used_tramadol)
269/(1231+269) #17.9 percent responded yes to: have you or anyone you know 
#               personally (friend or family) used Tramadol as a recreational drug?

# Cannibis Usage
table(egy_qual$cannabis)
egy_qual$cannabis[egy_qual$cannabis == ""] <- NA
egy_qual$used_cannabis <- ifelse(egy_qual$cannabis=="نعم",1,0)
table(egy_qual$used_cannabis)
287/(1213+287) #19.1 percent responded yes to: have you or anyone you know 
#               personally (friend or family) used cannabis as a recreational drug?

# Any Drug Usage

egy_qual$any_drug <- ifelse(egy_qual$used_cannabis==1,1,
                            ifelse(egy_qual$used_tramadol==1,1,0))
table(egy_qual$any_drug)
392/1500

# table
table(egy_qual$used_cannabis,egy_qual$used_tramadol)

egy_qual_na <- egy_qual %>% dplyr::select(starts_with("used")) %>%
  na.omit()

drug_table <- table(egy_qual_na$used_cannabis,egy_qual_na$used_tramadol)
drug_pro <- prop.table(drug_table)
drug_pro
#  73.9 have not or do not know someone for either
#   8.2 have tried/know someone who tried cannabis but not tramadol
#   7.0 have tried/know someone who tried tramadol but not cannabis
#  10.9 have tried/know someone who tried both
#  26.1 have tried/know someone who tried any or both



#########################################################################   
#                  DV: Drug Punishment                                  #   
#########################################################################

table(egy_qual$Drug_Punishment)
egy_qual$Drug_Punishment[egy_qual$Drug_Punishment == ""] <- NA
egy_qual$drug_punish_always <- ifelse(egy_qual$Drug_Punishment=="يجب أن يقضي الأشخاص الذين يتم إلقاء القبض عليهم بسبب حيازتهم للقنب الهندي (”حشيش“) عقوبة بالسجن الإلزامي",1,0)
table(egy_qual$drug_punish_always)

egy_qual$drug_punish_judges <- ifelse(egy_qual$Drug_Punishment=="يجب أن يحدد القضاة أحكام السجن على أساس كل حالة على حدة" ,1,0)
table(egy_qual$drug_punish_judges)

egy_qual$drug_punish_first <- ifelse(egy_qual$Drug_Punishment=="الأشخاص الذين اعتقلواوا على إثر حيازتهم للقنب الهندي (”حشيش“) لا ينبغي أن يقضوا عقوبة بالسجن الإزامي على الأقل في المرة الأولى والثانية\n",1,0)
table(egy_qual$drug_punish_first)

egy_qual$drug_punish_decrim <- ifelse(egy_qual$Drug_Punishment %in% 
                                        grep("يجب أن لا تعتبر حيازة القنب",
                                             egy_qual$Drug_Punishment, value=TRUE), 1, 
                                      ifelse(is.na(egy_qual$Drug_Punishment), NA, 0))
table(egy_qual$drug_punish_decrim)

egy_qual$Drug_Punishment_Policy <- ifelse(egy_qual$drug_punish_decrim==1,1,
                                          ifelse(egy_qual$drug_punish_first==1,2,
                                                 ifelse(egy_qual$drug_punish_judges==1,3,
                                                        ifelse(egy_qual$drug_punish_always ==1,4,NA))))
table(egy_qual$Drug_Punishment_Policy)


egy_qual$Drug_Punishment_Policy2 <- NA
egy_qual$Drug_Punishment_Policy2[egy_qual$Drug_Punishment_Policy==1] <- "Decriminalization of cannabis possession"
egy_qual$Drug_Punishment_Policy2[egy_qual$Drug_Punishment_Policy==2] <- "First two offenses have no jail time"
egy_qual$Drug_Punishment_Policy2[egy_qual$Drug_Punishment_Policy==3] <- "Judges evaluate on case-by-case basis"
egy_qual$Drug_Punishment_Policy2[egy_qual$Drug_Punishment_Policy==4] <- "Always mandatory prison sentences"

egy_qual$Drug_Punishment_Policy2 <- factor(egy_qual$Drug_Punishment_Policy2, levels = c("Decriminalization of cannabis possession",
                                                                                        "First two offenses have no jail time",
                                                                                        "Judges evaluate on case-by-case basis",
                                                                                        "Always mandatory prison sentences"))
table(egy_qual$Drug_Punishment_Policy2)

egy_qual$drug_depenalization <-  NA
egy_qual$drug_depenalization[egy_qual$Drug_Punishment_Policy==1] <- 1
egy_qual$drug_depenalization[egy_qual$Drug_Punishment_Policy==2] <- 1
egy_qual$drug_depenalization[egy_qual$Drug_Punishment_Policy==3] <- 0
egy_qual$drug_depenalization[egy_qual$Drug_Punishment_Policy==4] <- 0
table(egy_qual$drug_depenalization)


#########################################################################   
#                  Egypt: Drug Punishment                               #   
#########################################################################

egy_qual_punishment <- egy_qual[is.na(egy_qual$Drug_Punishment_Policy2)==F,]

drug_punishment <- ggplot(egy_qual_punishment, aes(Drug_Punishment_Policy2), na.rm = T) + 
  coord_flip() +
  geom_bar(aes(y = (..count..)/1486)) + 
  labs(title = "Policy Preferences for Drug Possession Punishment", 
       subtitle = "Egypt, March 2020 (Qualtrics Sample)", 
       y = "Proportion of Respondents", x = "Policy") +
  scale_y_continuous(expand = c(0,0)) +
  scale_x_discrete(labels = c("Decriminalization of possession", "No jail for first two offenses", 
                              "Case-by-case basis","Mandatory sentences")) +
  theme_bw() +
  theme(axis.text.x=element_text(size=9),
        axis.title.x=element_text(size=11),
        axis.text.y=element_text(size=9),
        axis.title.y=element_text(size=11),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())

drug_punishment


#########################################################################   
#                  Egypt: Drug Punishment & Drug Use/Exposure           #   
#########################################################################

drug_use_ex <- table(egy_qual_punishment$any_drug, egy_qual_punishment$Drug_Punishment_Policy2) 
drug_use_pro <- prop.table(drug_use_ex, 1) 
drug_use_pro <- as.data.frame(drug_use_pro)
colnames(drug_use_pro) <- c("drug_use_pro", "answer","propor")
drug_use_pro$percent <- (drug_use_pro$propor)*100


drug_use_likert_punish <- ggplot(drug_use_pro, aes(fill=answer, y=percent, x=drug_use_pro)) + 
  geom_bar(position="stack", stat="identity") +
  coord_flip() +
  labs(title = "Drug Use Penalization Attitudes", 
       subtitle = "Which is closest to your opinion regarding punishment for cannabis possession?", 
       y = "Percent of Sample", x = "", 
       fill = "") +
  scale_y_continuous(limits=c(-2, 102),
                     expand = c(0,0), 
                     breaks = c(0,25,50,75,100), 
                     labels = c("0%","25%","50%","75%","100%")) +
  scale_x_discrete(labels = c("No drug use or\nexposure to use","Used drugs or knows\nothers who have")) +
  scale_fill_grey(start=0.8, end=0.2,guide = guide_legend(reverse = TRUE)) +
  theme_bw()+theme(legend.position = "bottom", legend.direction = "vertical")

drug_use_likert_punish

png(file = "~/Dropbox/Egypt_Alex Sarah Liz/Figures/drug_use_likert_punish.png",
    width = 9, height = 6, units = 'in', res = 300)
drug_use_likert_punish
dev.off()


#########################################################################   
#                  DV: Drug Treatment                                   #   
#########################################################################


table(egy_qual$Drug_Treatment)
egy_qual$Drug_Treatment[egy_qual$Drug_Treatment == ""] <- NA

egy_qual$drug_treatment_decrease <- ifelse(egy_qual$Drug_Treatment=="لا، يجب على الحكومة المصرية تقليص دعمها المالي لمرافق معالجة الإدمان.",1,0)
table(egy_qual$drug_treatment_decrease)

egy_qual$drug_treatment_increase <- ifelse(egy_qual$Drug_Treatment=="نعم، يجب على الحكومة المصرية تقديم المزيد من الدعم المالي لمرافق معالجة الإدمان.",1,0)
table(egy_qual$drug_treatment_increase)


egy_qual$drug_treatment_same <- ifelse(egy_qual$Drug_Treatment=="يجب أن يظل الدعم الحكومي لمرافق معالجة الإدمان على حاله دون تغيير.",1,0)
table(egy_qual$drug_treatment_same)


egy_qual$Drug_Treatment_Policy <- ifelse(egy_qual$drug_treatment_decrease==1,1,
                                         ifelse(egy_qual$drug_treatment_same==1,2,
                                                ifelse(egy_qual$drug_treatment_increase==1,3,NA)))

table(egy_qual$Drug_Treatment_Policy)

egy_qual$Drug_Treatment_Policy2 <- NA
egy_qual$Drug_Treatment_Policy2[egy_qual$Drug_Treatment_Policy==1] <- "Decrease funding"
egy_qual$Drug_Treatment_Policy2[egy_qual$Drug_Treatment_Policy==2] <- "Keep funding the same"
egy_qual$Drug_Treatment_Policy2[egy_qual$Drug_Treatment_Policy==3] <- "Increase funding"

egy_qual$Drug_Treatment_Policy2 <- factor(egy_qual$Drug_Treatment_Policy2, levels = c("Decrease funding",
                                                                                      "Keep funding the same",
                                                                                      "Increase funding"))

table(egy_qual$Drug_Treatment_Policy2)


#########################################################################   
#                  Egypt: Drug Treatment                                #   
#########################################################################

egy_qual_treatment <- egy_qual[is.na(egy_qual$Drug_Treatment_Policy2)==F,]

drug_treatment <- ggplot(egy_qual_treatment, aes(Drug_Treatment_Policy2), na.rm = T) + 
  coord_flip() +
  geom_bar(aes(y = (..count..)/1476)) + 
  labs(title = "Policy Preferences for Drug Treatment Funding", 
       subtitle = "Egypt, March 2020 (Qualtrics Sample)", 
       y = "Proportion of Respondents", x = "Policy") +
  scale_y_continuous(expand = c(0,0)) +
  scale_x_discrete(labels = c("Decrease funding", "Keep funding the same", 
                              "Increase funding")) +
  theme_bw() +
  theme(axis.text.x=element_text(size=9),
        axis.title.x=element_text(size=11),
        axis.text.y=element_text(size=9),
        axis.title.y=element_text(size=11),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())


drug_treatment


#########################################################################   
#                  Egypt: Drug Treatment & Drug Use/Exposure            #   
#########################################################################

drug_use_ex <- table(egy_qual_treatment$any_drug, egy_qual_treatment$Drug_Treatment_Policy2) 
drug_use_pro <- prop.table(drug_use_ex, 1) 
drug_use_pro <- as.data.frame(drug_use_pro)
colnames(drug_use_pro) <- c("drug_use_pro", "answer","propor")
drug_use_pro$percent <- (drug_use_pro$propor)*100


drug_use_likert_treat <- ggplot(drug_use_pro, aes(fill=answer, y=percent, x=drug_use_pro)) + 
  geom_bar(position="stack", stat="identity") +
  coord_flip() +
  labs(title = "Drug Treatment Funding Attitudes", 
       subtitle = "Which is closest to your opinion regarding government funding for drug treatment?", 
       y = "Percent of Sample", x = "", 
       fill = "") +
  scale_y_continuous(limits=c(-2, 102),
                     expand = c(0,0), 
                     breaks = c(0,25,50,75,100), 
                     labels = c("0%","25%","50%","75%","100%")) +
  scale_x_discrete(labels = c("No drug use or\nexposure to use","Used drugs or knows\nothers who have")) +
  scale_fill_grey(start=0.8, end=0.2,guide = guide_legend(reverse = TRUE)) +
  theme_bw()+theme(legend.position = "bottom", legend.direction = "vertical")

drug_use_likert_treat

png(file = "~/Dropbox/Egypt_Alex Sarah Liz/Figures/drug_use_likert_treat.png",
    width = 9, height = 6, units = 'in', res = 300)
drug_use_likert_treat
dev.off()



###############################################################################
#                      ANALYSIS FOR APPENDIX                                  #
#                                                                             #
###############################################################################

# Create list of model names
models.qual <- list()

# Run models
models.qual[[1]] <- lm(drug_depenalization~any_drug + resp_fem + age + muslim, data=egy_qual)

models.qual[[2]] <- lm(drug_depenalization~any_drug +  employed + unemployed + secondary_plus +
                         resp_fem + age + muslim, data=egy_qual)

models.qual[[3]] <- lm(drug_treatment_increase~any_drug + resp_fem + age + muslim, data=egy_qual)

models.qual[[4]] <- lm(drug_treatment_increase~any_drug +  employed + unemployed + secondary_plus +
                         resp_fem + age + muslim, data=egy_qual)

names(models.qual) <- c("Drug Depenalization", "Drug Depenalization", "Drug Treatment", "Drug Treatment")

stargazer(models.qual, 
          out = "Tables/drug_qual.tex",
          title = "Qualtrics Survey: Exposure to Drug Use and Policy Attitudes",
          label = "drug_qual",
          type = "text", 
          omit.stat = c("f","adj.rsq","ser"), 
          #omit = c("ind_ID"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.qual),
          dep.var.labels = c("Probability of Support for Policy"),
          #covariate.labels = c(), 
          add.lines = list(c("Sex/Age/Religion Controls", rep("Y", length(models.qual))),
                           c("Additional Controls", "N","Y","N","Y")), 
          notes.align = "l",
          notes = c("OLS models. All models control for respondent sex, age, and religion. ",
                    "Models 2 and 4 also control for employment and education."))



